RADIATION PRESSURE SUPPORTED AGN TORI WITH 
HARD X-RAY AND STELLAR HEATING 
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ABSTRACT 

The dynamics and structure of toroidal obscuration around AGN remain 
uncertain and controversial. In this paper we extend earlier work on the dy- 
namical role of infrared radiation pressure by adding the effects of two kinds 
of distributed heating: Compton-heating due to hard X-rays from the nucleus 
and local starlight heating. We find numerical solutions to the axisymmetric hy- 
drostatic equilibrium, energy balance, and photon diffusion equations including 
these effects. Within the regime of typical parameters, the two different sources 
of additional heating have very similar effects: the density profile within the 
torus becomes shallower both radially and vertically, but for plausible heating 
rates, there is only minor change (relative to the source-free case) in the distri- 
bution of column density with solid angle. The most interesting consequence of 
distributed heating is that it selects out a relatively narrow range of parameters 
permitting an equilibrium, particularly (L/Le)/tt- We discuss the implications 
of both the narrowness of the permitted range and its approximate coincidence 
with the range inferred from observations. 

Subject headings: 



1. INTRODUCTION 



The obscuring torus is one of the key components to the aniso tropic appearanc e of 
AGN. Althoug h much observational evidence exists to directly (e.g. IJaffe et al.l l2004f) or 
indirectly (e.g. lBarthellll989l ; Idi Serego Alighieri et.al.l Il994l ; IZakamska et al.l 120051 ) confirm 
the existence of this structure, there is little understanding of its dynamics. The central 
question is the nature of the mechanism that supports the torus's large geometrical thickness 
against gravity. 

Numerous ideas have been proposed to answ er this question. The first suggestion was the 
clumpy torus model (IKrolik fc Begelmanlll988l ). In that model, the torus consists of highly 
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clumped gas and dust, and the clumps undergo highly supersonic motions. To avoid rapid 
collisional loss of kine tic energy, a large magn etic field is needed to ensure that each collision 
is sufficiently elastic. iKonigl fc Kartjd (119941 ) presented another possible approach. In their 
model, they argued that a magneto-centrifugal wind could replace the torus. However, this 
model faces difficulties to explain the origin of the large-scale magnetic field and the source 
of the large energy need ed to drive the wind . Another possibility is that the support is 
from radiation pressure (IPier &: Krolikl Il992al ) . The optical through soft X-ray continuum 
of the nucleus is absorbed and re-emitted in infrared light by dust at the inner edge of the 
torus; then the large opacity in that band couples the radiation flux to the ga s and p r ovides 
a strong radiation force to balance the gravity. Following this basic idea, iKrolikl (120071 . 
hereafter K07) constructed an idealized model, and via this model found self-consistent 
hydrostatic equilibrium solutions analytically. These solutions demonstrate that infrared 
radiation pressure is able to support the geometrically thick structures around AGN. For 
simplicity, that work did not consider any internal sources of heating, such as the Compton 
heating due to hard X-rays penetrating the torus interior, or the heating from local starlight 
irradiating the dust. Both of them contribute a positive divergence of flux to the energy 
equation, which can strongly affect the configuration of the torus and even the existence of 
equilibrium solutions. 

It is the object of this paper to construct a generalized radiation support model by 
including these local heating mechanisms. We first construct the physical model in section 
[21 introducing the basic equations and assumptions adopted in this work. Section [3] shows 
how we solve these equations. After defining dimensionless parameters and identifying ap- 
propriate boundary conditions, we describe in detail the numerical method implemented in 
this work. The results and discussion are presented in section [H and the conclusions follow 
in section [51 



2. THE PHYSICAL MODEL 

We choose 2-d axisymmetric geometry to explore this picture. All physical quantities 
are written in cylindrical coordinates on the r — z plane. To be appropriate to flattened 
geometries, three simplifying assumptions are adopted. First, we take Q as the local or- 
bital frequency, which at all heights z equals the rotation rate of a circular orbit in the 
torus midplane at radius r. Second, we follow only the component of angular momentum 
parallel to the torus axis, and we assume that the gas's specific angular momentum has 
magnitude jr 2 Q with j = j(r, z) < 1. Third, the radial and vertical components of gravity 
are approximated by rVL 2 and zQ 2 . In fact, Q(r, \z\ > 0) < Q(r, z = 0), so this approxima- 
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tion slightly overestimates the strength of gravity. For example, in a point mass potential, 
2 (r, z) = GM BH /(r 2 + z 2 ) 3//2 < f2 2 (r, z = 0), where M BH denotes the mass of the central 
black hole. We also assume the interior of the torus is in hydrostatic equilibrium, so that: 

k¥/c = -g e// = rfi 2 (l - j 2 )e r + zQ 2 e z , (1) 

where the infrared radiation flux is F, k is the opacity per unit mass, and g e ff is the net 
gravity. 

Instead of solving a complete transfer problem at all relevant frequencies for all pho- 
ton directions, we approximate the radiation flux by solving the diffusion equation with a 
thermally-averaged opacity. In this approximation, the flux is obtained from the gradient of 
the radiation energy density: 

F = -^-VE, (2) 
where p is the gas mass density and E is the radiation energy density. 

If the only source of infrared radiation is the conversion via dust reradiation of optical 
and UV photons at the inner edge of the torus, then in the body of the torus 

V ■ F = 0. (3) 

However, the existence of distributed sources in the torus is also possible. For instance, when 
hard X-rays penetrate d eeply into the torus material, local heating due to Compton recoil 
(e.g., IChang et al.l 120071 ) can be considerable. It is also possible that local star formation 



is sufficiently strong; th at stellar luminosity may supplement the AGN's radiation force (see 



Thompson et al.ll2005l ). 



In this paper, we explore both of these. In the former case, if the torus is optically thin 
to hard X-rays, a more general formula to describe the energy conservation reads 

v - F = ^WT^)"' aTf " (4) 

where Lx is the luminosity in hard X-rays, n e is the electron number density, and f c is the 
ratio of the energy gained by electrons during each collision to the photon energy. Even 
bound electrons behave as if they are free when scattering X-rays with energy greater than 



roughly 3-4 keV (jKrolik!ll999l ). and Klein-Nishina effects are negligibly small for hard X-ray 



photons < lOOkeV, thus ax, the Thomson cross section, is the appropriate cross section. 

We have much more freedom to choose the distribution of internal starlight. A reason- 
able assumption is to adopt the Schmidt Law that the star formation rate is proportional to 
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the 3/2-power of the gas density jKennicuttl[l998h . and assume that the local stellar lumi- 
nosity is proportional to the star formation rate. Due to the large optical depth to the UV 
and optical, the stellar radiation is assumed to be absorbed in situ and reproduced in the 
infrared. As in equation HI we can write the equation of energy conservation including local 
stellar heating as 

/ \ 3 / 2 

(5) 



V-F 



Pin 



where S is a coefficient with units of erg cm 3 s 1 which describes the strength of the sources, 



and p in = p(r in , 0) is the gas density along the inner edge of the torus (r 
the midplane. 



measured on 



With equations fT|2l and either H] or [5j it is possible to solve for all three unknowns, F, 
E and p. 



THE GENERALIZED SOLUTIONS 



3.1. Preparatory Work 

To find the energy density and matter density from the three equations introduced in 
the last section, we need to combine them and simplify them. 



Putting the flux equation [2] together with the hydrostatic equilibrium equation [TJ we 



have 



3p 



VE = rQ 2 (l - f)e r + zQ 2 e z 



(6) 



which relates the energy density to the dynamics. 



Combining equation [T] and either [4] or [5], one gets the relation between the dynamics 
and the local sources of heat: 



V- |^[rfi 2 (l -j 2 )e r + ^ 2 e,]} = R, 



(7) 



where R = n e aTf c Lx/^(r 2 + z 2 ) in the case of hard X-ray heating and R = S (p/Pi n ) 3 ^ 2 in 
the st ellar heating case. According to the most recent dust opacity models (e.g. JSemenov et al 
20031 ). the Rosseland mean opacity is a mildly changing function of the temperature in the 
range 100 — 1000 K, which is also the interior temperature of the obscuring tori as found 
by detailed radiation transfer studies (IPier &: Krolikl Il992bl ; lEfstathiou &: Rowan- Robinson 



1995: Granato et al. 1997: Nenkova et al 



20021 ). On this ground, we approximate k as con- 
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stant. Equation [7] then reduces to 

r% + 2(1 - <*)j 2 + (2a - 3) = -^R, (8) 

where a is defined by O(r) = f2j n (r/rj n ) _a in order to allow for more general potentials than 
that of a simple point-mass. VL in = Q(r in ) in that definition represents the orbital frequency 
measured on the inner edge of the torus. Without terms on the right hand side, we can solve 
equation M easily because there is no longer any dependence on z. When R ^ 0, we can 
treat the right hand side as a perturbation and solve the exact equation iteratively. Detailed 
description of this method will be given in subsection 13.41 

Once j 2 (r, z) is found, we can turn back to equation [6] and separate it into two equations: 

1 dE 1 dE 

P ~ ~3zQ 2 ~dz ~ ~3rQ 2 {l-j 2 )~dr~' ^ ' 

The second equality in the above equation allows us to rewrite the partial differential equation 
in characteristic form 

dE _ dE dz + dE dr _ Q 
ds dz ds dr ds 

with 

dz 1 dr 1 
ds z ' ds r(l — j 2 ) 
^From the characteristic form, we find E is constant along contours parameterized by A: 

\z 2 + f dr'r'[l-f(r',z)}=\, (12) 

where r* is arbitrary. 

However, to find the exact values of E = E(r, z) at distinct locations, we need to know 
the energy density E(X) along any path on which A varies monotonically. For convenience, 
we can pick the path to run outward along the r axis starting from the radius of the inner 
edge r in . We then have 

dE = (mdr_ = _ 3 ^ 2 ^ 
dX dr dX 

In order to achieve a solution, this equation requires advance knowledge of p(A) on its 
path. It is convenient in this context (in which we have already written Q oc r~ a ) to 
consider density boundary conditions that also have power-law dependence on radius, i.e. 
p(r, 0) = Pin(r/rj n )~ 7 . With this choice of gas density, the energy density can be easily found 
by integrating 

rlR(r (X\ 

= -3p(r,0)rn 2 [l-f(r,z)}, (14) 



dr 
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that is 

E{r, 0) = E in - I 3p(r', 0)r'ft 2 [l - j 2 (r', 0)}dr', (15) 

where E in is the energy density measured at point (r in ,0). Thus, given j 2 , which is the 
solution of equation [HJ the density along the equatorial plane, and the energy density on the 
equatorial plane, the complete solution for E can be achieved with equations [12] and [15j The 
matter density p = p(r, z) can then be computed from either equality of equation 



3.2. Governing Parameters 

We have mentioned two important dimensionless parameters in the previous subsection: 
a and 7. They determine the shape of the gravitational potential and the density profile 
on the equatorial plane. Besides those two, we still need several others to parameterize our 
problem. One of these is j in , defined as j in = j(r in , 0), which indicates the rotational support 
at the inner edge. Greater j in means that the torus requires a smaller radial thickness to 
reach the full Keplerian angular momentum. Another parameter is r* = K,pi n ri n , which sets 
the optical depth scale. If the density declines outward, r* must be at least several to satisfy 
the diffusion approximation. Because the Rosseland mean opacity of dust per unit mass of 
gas is ~ 10 - 30 times as great as the Tho mson opacity per unit mass for temperature in 



the range 100-1000 K (ISemenov et al.ll2003l ). tt = ^TPinfin ~ (0.03 — 0.1)r*. For simplicity, 
tt = 0.05r* and r* = 10 (so that tt = 0.5) are adopted in this paper if we do not say 
otherwise. 

We also need a parameter Q = 3p in r 2 n fl 2 n / E in to relate the orbital energy to the radia- 
tion energy. If no local heating exists in the torus, or the extra sources are negligibly weak, 
then the only contribution to Ei n is the absorbed radiation of the nucleus, which is mainly 
in the UV band. Thus we can write E in ~ L uv h/(4:iir 2 n c), where L uv denotes the luminosity 
of the nucleus in the UV, and h is a blanketing factor, telling us by what factor the inner 
edge energy density is enhanced compared to what it would be in vacuum. However, the 
local heating also contributes to the energy density, and therefore we have 

L uv h ff rRe~ T ( r ' z ^ L vv h 
E in = -j— 2- + / / drdz- = j iL 2-/*- (16) 

The integral over the torus gives the contribution to the energy density due to the local 
sources. Here we assume the gas is axisymmetrically distributed and r(r, z) is the infrared 
optical depth from the location of the source to the position (r in , 0). /* is a correction factor 
> 1. With Ei n defined this way, the parameter Q can be easily rewritten in terms of more 
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familiar quantities: 



Q 



r* M(< r in ) k t L E 



(17) 



where M(< rj n ) is the total mass interior to rj n , and Le is the Eddington luminosity. As 
discussed in K07, if r* ~ 10-30, and the blanketing factor h > 2, then reasonable values of 
Q would be ~ 0.1 - 10. 

The last two parameters characterize the internal sources. Firstly, the X-ray heating 
requires a parameter X = Lxf c /{Luvf*)- Combining this definition with equation [T6l we 
find that the correction factor /* is 



T T X 



where/; = ff drdz(p/ p in )(r/r in )e r{r ' z) /{[(r-r in ) 2 + z 2 ][(r/r in ) 2 + (z/r in ) 2 }}. Consequently, 
in terms of X, the ratio of the hard X-ray luminosity to the total luminosity from the nucleus 
is: 

(19) 



J X 



L X + L L 



In AGN hard X-ray sp ectra, the photon spectral index ranges from ~ 1-3 (IBeckmann et al. 
20061 : iTozzi et al.l 120061 ) ; the averaged fractional energy lost f c =<hv /m P c 2 > is then ~ 0.1. 
Using the fact that bolometric corrections for 2 — 10 keV X-rays are ~ 8-60 (IMarconi et al. 
20041 ). assuming a photon index of 2, and extrapolating the hard X-ray spectrum up to 100 
keV, we estimate that the ratio of total X-ray luminosity to bolometric Lx/L ~ 0.1-0.3. We 
therefore expect X < 0.05. 

Secondly, the dimensionless parameter that determines the internal stellar heating is 
defined as P = 127rr J 3 n S'/(L (7V / ;)! ). After substituting for S in equation [TBI with this definition 
of P, we obtain the function /* for stellar heating: 



P 

1 vh 

3h 2 



-i 



(20) 



where h — Jj drdz (p/pi n ) 3 ^ 2 (r/r in )e r ( r > 2 )/[(r — r in ) 2 + z 2 ]. With /* known, we can relate 



S to P, and then calculate the luminosity of the starlight L star 
Its ratio to the AGN luminosity is then 



■'star 



UP 



3/2 



JJdrdzA7rrS(p/p m ) 3 / 2 . 



(21) 



Note that L star cannot be observed directly due to the large infrared o ptical depth in the 
torus. Recent work based on integral field spectroscopy with SINFONI ( jDavies et al.l 120071 ) 
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provides a conservative estimate for the stellar luminosity within the central 10 pc: the ratio 
of starlight to AGN light is less than a few percent. Therefore, the possible values of L star /L 
could be ~ 0.001 — 0.01, and from this ratio we can constrain the values of P. 

With all of these typical scalings and parameters, we can now put the principal equations 
(eq. [BJ [T2land [T5l) in non-dimensional form. For example, equation [8] can be rewritten as 



r 

Tin 



df 



where 



for X-ray heating, and 



r 



d(r/r, 
3r*r T 



+ 2(1 - a)f + (2a - 3) = T, 



XI — 



Qh 



r = -— -p 



{r/r in ) 



2a 



Qh' 



_P_ 

Pin 



_{r/r in ) 2 + (z/r in ) 2 

3/2 



r 

? in 



2o 



(22) 
(23) 

(24) 



for stellar heating. Because the parameter h appears in the perturbations only in combination 
with Q, it is convenient to absorb the effect of h into Q. From here on out, we fix h = 5. 

In sum, we have six parameters that govern the character of the solution: j in , a, 7, Q, r*, 
and X or P. All are independent except 7. Since it goes into the density boundary condition 
along the equatorial plane, we discuss it in the next subsection on boundary conditions. 



3.3. Boundary Conditions 

Three boundaries were discussed in K07, and they also apply here. The first one is the 
inner boundary or the inner edge. Since we are interested in solutions in the interior of the 
torus, we require r > r in and simply choose a vertical inner edge r e d ge (z) = r in in this paper. 
As discussed at greater length in K07, the inner edge is not a physical boundary. All it does 
is mark the limit of the region within which we evaluate our solution; there is no reason 
to think that the actual inner edge of the cool, dusty material is exactly vertical. At this 
stage of our understanding, we choose to leave its actual shape unspecified for two reasons: 
One is that, by seeking hydrostatic solutions inside the torus, our problem is sufficiently 
determined mathematically as to obviate the need for another boundary condition on the 
inner edge. In K07, we showed that, to the degree that we can estimate the solution to the full 
radiation transfer problem inside the torus "hole", our solutions are at least approximately 
self-consistent. The other reason is that the physics determining the real shape of the inner 
edge is a complicated brew of hydrodynamics, photoionization physics, and dust-sputtering 
dynamics far beyond the scope of this simplified model. 
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Secondly, the energy density should not go negative, which places a constraint on the 
outer boundary of the torus. Anywhere E < 0, the solution is unphysical and must therefore 
be discarded. Thirdly, the photosphere acts as a much stricter outer boundary because 
the diffusion approximation is valid only in the optically thick region. At the photosphere, 
where t z = J 2 °° dz' np(r, z') = 1, the estimated diffusive flux should roughly match the flux 
as evaluated in the free-streaming case: 

~cE(r,z). (25) 

As discussed in K07, we expect this boundary condition to be satisfied only to within a 
factor of 3. 

Another factor that plays a role in determining the outer boundary is the requirement 
that j(r, z) < 1: greater j would make hydrostatic equilibrium impossible. Although there is 
no physical inconsistency in positing a sub-Keplerian outer boundary, it is hard to understand 
the dynamical state of the matter beyond this edge. What is the supporting force outside 
that boundary? Suppose that it is rotationally supported, what then makes the transition 
from partial radiation force support to full rotational support? For this reason, we define 
the outer boundary r max by requiring j(r max , 0) = 1. 

In practise, given a set of parameters of which all except 7 are fixed, the requirement 
that E > coupled with the photospheric boundary condition determines the proper value 
of 7. In this sense, 7 is a sort of eigenvalue, and the density profile on the midplane p(r, 0) = 
Pin{ r / r in)~ 1 is not an independent boundary condition. We also find that the photospheric 
boundary condition is best matched at the smallest 7 such that E(r, 0) > everywhere in 
the range r in < r < r max . 



3.4. Numerical Approach 

The basic equations listed in section [37X1 cannot be solved analytically, so we must invent 
a numerical method. Let us begin by considering equation [22J The difficulty in solving it 
comes from the inhomogeneous terms, which depend on p(r, z), a quantity we know only 
after solving the problem. If the right hand side is zero, however, the equation reduces to 
an ordinary differential equation in j 2 for which a solution can be easily found. Right hand 
sides that are "small" can therefore be regarded as perturbations. Once the zeroth order 
solution (the homogeneous solution) has been obtained, we can find the first order solution 
by restoring the right hand side and using p(r, z) from the zeroth order solution as the initial 
guess. Following this procedure and iterating, the rt-th order solution can be solved if the 
(n — l)-th order solution is in hand. Because we do expect the internal heating effects to be, 
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in some sense, "small" , we may reasonably hope for convergence. 

K07 has already given a detailed description of how to solve the homogeneous version 
of equation [22] Here we would like to review this method briefly. After neglecting the right 
hand side of equation [22], and treating j 2 as a function of r only (because there is no longer 
any z dependence left), its solution is: 

i( 2 o)(r) = til + f(a)](r/r in r^ _ f{a) (26 ) 

for a 7^ 1, where f(a) = 0.5(3 — 2a) /(a — 1), and the subscript in parenthesis denotes the 
order of the solution. For a = 1, j? -\(r) = jf n + ln(r/r in ). 

With j? Q s (r) known, we can substitute it into the characteristic curves of energy density 
E in equation [12], giving 

K^) 2+ i<^T)fc) 2 -^ +/(a,, fc)^ A (27) 

if a ^ 1 and 



1 / z 



2-1 +;(!-£,)(- 



1 / r 

2 V Tim 



In I — 



A 



(2? 



for a = 1. 

Similarly, we can easily integrate equation [TH and write out the integral in equation [15] 



as 



£ (0) (r,0) = EiJl-Q 



2 -2a - 7 



2-2a-7 



+ 



j1 + /(a) 

7 



for a 7^ 1 and 



s (0) (r,o) = ^Ji-g 



r 



-7 



in / 7 
1 



— - lnr-l+^+- 



7 



7 V 7 



if a; = 1. 



(29) 



(30) 



Finally the zeroth order energy density throughout the plane can be found by using 
the explicit form of Ery on the midplane and the equivalence contour of Er ) in r — z space 
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from equations [27] or [28] Then the mass density p can be computed by either equality in 
equation [9j 

The next step is to compute the first order solution. Notice that now j 2 in equation I2"2l is 
a function of both r and z. The local heating actually changes the distribution of the angular 
momentum, and the contours of j 2 are no longer vertical lines, but are instead shifted and 
bent away from the inner edge near the equatorial plane. Integrating equation [22} we have 



j(i)(r, z) = J ( 2 0) (r, z) + |y r (0) (r', z)r' x - 2 «dr> j r 2 ^\ (31) 

where T( ) (r, z) is the right hand side in that equation, which is also the perturbation evalu- 
ated with the zeroth order solution of p. Remember that the parameter Q in the perturbation 
is the same as that in the zeroth solution. The next step after we obtain the angular mo- 
mentum distribution is to recalculate the characteristic curves of the energy density. By 
performing the integration in equation [12], we can numerically find A(r, z). 

The values of E and its corresponding characteristic parameter A at the boundary are 
required in order to visualize the contours of the perturbed energy density. Equation [15] 
determines the equatorial values of E: 



E (1) (r, 0) = E in \ 1 - Q I ( — ) 7 [l - jf^r', 0)] — } , (32) 



where we have fed in the boundary condition for mass density, i.e. p(r, 0) = ft„(r/r ffl )' 7 . As 
discussed in subsection 13.3] the requirements jm(r = r max , 0) = 1 and E(X)(r < r max , 0) > 
help us find the proper 7. 

Finally, we obtain the first order solution of the radiation energy density in the torus 
by interpolating on the r — z plane. Again the mass density is determined by the partial 
differential equation O 

To achieve a higher order of accuracy, we can follow the whole procedure again by 
substituting the lower order solutions into the complete set of partial differential equations 
and keep iterating. We terminate the procedure when further iterations no longer change the 
solution. Lastly, after the iterations have converged, we test whether the solution satisfies 
the photospheric boundary condition, accepting the result only if it does. A flow chart (Fig. 
[Q is presented to illustrate the procedure more clearly. 

For all the numerical calculations, we construct an evenly spaced grid to cover the 
region r in < r < r max and < z < z max , where r max and z max are determined by the 



boundary conditions. We use the sum Yl 



< N(t)E as the convergence criterion, 
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where the sum includes only those grid points at which Eha > 0, and NrjA is the number of 
those points. In addition to that criterion, we insist that the sum over those points should 
decrease monotonically as the order of the solution k grows. We set e = 10~ 2 , but much 
higher relative accuracies (often 10~ 4 ) can be reached after several steps of iteration. The 
fact that the iterative method is strongly convergent proves the numerical approach based 
on the perturbative approximation is valid. 

4. RESULTS AND DISCUSSION 

4.1. Examples of Typical Solutions 

Let's find the solution for typical parameters ji n = 0.5, a = 1.5, r* = 10, Q = 4 and 
X = 0.02. Those parameters describe a torus half-supported rotationally at the inner edge, 
in a point mass potential, with a column density ~ 10 24 cm" 2 in the midplane, L ~ O.lLEdd, 
hard X-ray luminosity ~ 16% of the total, and f c = 0.1. The 7 determined by best fitting the 
boundary conditions is ~ 1.43. A similar case with X replaced by P = 0.025 also requires 
7 ~ 1.43; the corresponding stellar luminosity is ~ 6 x 10" 3 L. These two solutions are 
presented in Figures [2] and [3J Like the unperturbed (no local sources) solutions shown in 
K07, the contours of the radiation energy density for the generalized solutions are extended 
upward, and the contours of the constant density are extended radially. The white curve in 
both figures shows the photosphere. Our diffusion approximation is validated by the fact 
that most of mass of the torus is within the optically thick region. 

In both cases, the correction factor /* is almost unity, so both solutions have the same 
L uv . We compare them in detail in Figure HI The distributions of energy density and matter 
density are nearly the same in the optically thick zone, which suggests a rough mapping 
between these two local heating cases. In other words, if a solution with one internal heat 
source is found within the proper parameter space, a very similar solution with the other 
must also exist. The slight distinction at larger distance is due to the different radial and 
vertical dependence of their perturbations in equation [221 

4.2. Comparison With the Unperturbed 

To better illustrate the effects of the local sources, we consider larger X and P. Given 
j in = 0.5, a = 1.5, r* = 10, Q = 4, and X = 0.06 (L x /L ~ 0.38) or P = 0.05 (L star /L ~ 
0.014), we find that 7 = 1.5. Because of the similarity between the two internally-heated 
solutions, we need only compare one of them with the unperturbed. Here we choose X-ray 
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heating. The correction factor /* for X = 0.06 is ~ 1.05. Taking into account this factor, 
we find that an unperturbed solution with Q ~ 4.2 possesses the same L uv as that of the 
perturbed. 7 for this solution is 1.6. Because there is more support at large radius with X-ray 
heating, the density profile becomes flatter both radially and vertically than the source-free 
one (Fig. Ej). Meanwhile, X-ray heating also causes the energy density to decrease less 
rapidly away from (r in , 0) than in the unperturbed case because we are comparing at fixed 
L uv and there is now additional internally-generated flux due to the local heating. 

Further investigation of the distribution of j 2 provides a clearer picture of the perturbed 
and unperturbed solutions (Fig. [6]). In the interior of the torus near the inner edge, the 
infrared radiation pressure is large enough to balance gravity, so the presence of internal 
sources does not effect j 2 too much; however, at large radius, contributions from the local 
sources are relatively strong, while infrared flux from the inner edge diminishes. Particularly 
in the equatorial plane, the additional radiation support in the radial direction reduces the 
need for rotational support. As a consequence, the radial gradient of j 2 becomes shallower 
than in the case without local heating. 



4.3. Exploring Parameter Space 

Holding X or P fixed, the allowed solutions in the Q — 7 plane fall onto a track with 
small thickness. The thickness is due to the imprecision of the boundary condition required 
at the photosphere. Following the track, the parameter 7 grows as Q increases. There are 
no solutions above or below the track. Parameters in the region below it fail the outer 
boundary criterion that j = 1 at the maximum radius; those above the track do not satisfy 
the boundary condition at the photosphere. There is also a starting point for each track 
(Qmin,1min), such that no solutions can be found with smaller 7 and Q. This fact, too, is 
an example of converged solutions that fail the boundary condition on the photosphere. In 
particular, when Q < Q m i n , \F/cE\ is too small, which means gravity is too weak in the 
torus, so that no hydrostatic balance can be achieved. The starting point moves toward 
larger 7 and Q when X or P increases, while the track rises a bit due to the change of the 
energy density contributed from the local sources. This result is a corollary of the general 
picture we have presented: if UV-derived radiation support can, on its own, balance gravity, 
equilibrium in the presence of additional radiation force requires a smaller UV luminosity. 

As examples, we plot tracks with X = (P = 0), X = 0.01 (P = 5 x 10~ 3 ) and 
X = 0.10 (P = 0.07) in Figure [7^, which represent zero, weak and strong local sources 
respectively. We choose the parameters X and P so that the solution tracks with different 
local sources can be matched very well and only one curve is drawn for each set of X and 
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P. In K07, it was shown that 7 becomes unrealistically large for Q > 6, but Q m in could be 
as small as 0.1. In Figure [7^, we see that the perturbed solutions have a smaller range of Q. 
Even for X = 0.01 (equivalent to P = 0.005), the minimum Q permitting a solution is ~ 3. 
In the regime of typical parameters, this range corresponds to (L/Le)/tt ~ 0.1-0.2, with a 
tolerance of several. With X ^ 0.1 or P ^ 0.07, Q must be greater than 4. When we push 
the parameters to an unrealistic limit - X > 1.0 (P > 0.9), no solutions can be found with 
7 < 3. 

We also investigate how the introduction of non-zero X or P alters the solution's depen- 
dence on ji n and a. Like the starting points of the solution tracks in Q — 7 plane, the tracks 
in the ji n — 7 and a — 7 planes possess end points. With increasing luminosity of the internal 
sources, for fixed Q, the end points move toward smaller j in and larger 7 (Fig. [7b), and 
smaller a and larger 7 (Fig. [7b). That means, to find dynamical balance, a stronger source 
in the torus demands less rotational support at the inner edge, and a less steep gravitational 
potential profile in the interior. As a consequence, if a torus has large X or P, it must 
have a low orbital speed at rj n and a flat potential inside. The latter might be particularly 
compatible with large P, as the flattened potential presumably reflects the contribution of 
stellar mass. 

The last panel (d) in Figure [7J shows the roughly linear correlation between P and the 
luminosity ratio L star / L as a function of Q, whose proportionality coefficient is determined by 
the integral in equation [2B Because larger Q indicates less radiation support and therefore 
less matter density in the torus, the coefficient falls with increasing Q. Similar results can 
be found if the correlation is plotted as a function of ji n or a: the slopes are always positive, 
which means that a larger P indicates a larger stellar luminosity fraction. 

The characteristic optical depth r* enters in several distinct ways: with regard to the 
IR support due to converted UV radiation, the only effect it has independent of its presence 
in Q is to increase the opacity, so that the photosphere rises with increasing r* is all other 
parameters are held constant. Otherwise, an increase in r* is equivalent to an increase in 
Q. In the perturbations, however, it has a different effect: X-ray heating is proportional 
to t^t^X/Q, while Q oc r*, so that in one sense the heating rate is proportional to only a 
single power of the column density. On the other hand, if Q is held fixed, the heating rate 
is proportional to the square of the column density. Similarly, the stellar heating rate is oc 
r*P/Q, so that the combination t*/Q is independent of the column density, but this heating 
rate rises linearly with column density at fixed Q. When considering all of these scalings, it is 
important to note that both t? and T* are defined as characteristic optical depths {Kpi n r in ) 
rather than actual optical depths along any particular ray. For our assumption of X-ray 
free-streaming, the actual Thomson optical depth from the nucleus to any particular point 
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in the torus should be < 1. That requirement should always be satisfied if tt < 1; because 
the density falls off rapidly with increasing r and increasing \z\, it should be satisfied in the 
majority of the torus volume even when Tj- > 1. 



4.4. Comparison With Observations 

One measurable diagnostic of the torus is the column density of matter along the line 
of sight. Although it is difficult to measure the inclination angle of the torus to the line of 
sight in an individual object , one can still investigate the statistical distribution of column 



densities to obscured AGN (IRisaliti et al.lll999l ; iTriester et al.ll2004l ). This distribution can 



also be predicted by our model because the probability of a given column density is sim- 
ply proportional to the solid angle associated with the polar angle producing that column. 
However, there is a certain level of arbitrariness in this prediction due to the guessed shape 
of the inner edge. Nonetheless, if we assume a vertical inner edge, the generalized solutions 
predict marginally wider ranges of the column densities than the unperturbed due to extra 
radiation pressure support. Compared with the unperturbed, although most of the solid an- 
gle is still associated with the higher column densities, the shape of that distribution tends 
to be flatter. 

In Fi gure [HI we show the predicted distributions for different X and P , keeping L uv 
fixed. The descriptions of the curves and parameters are listed in Tabled] and [2j We find that 
in both cases the distribution gradually extends to higher column densities when stronger 
local sources are present; however, the change is very small and could be unmeasurable. 
For large enough X and P, the distribution curves roll over as they approach the largest 
column densities, an effect caused by the extended "foot" at large radius near the midplane 
in those solutions (e.g., the right panel of Fig. [3] and [5]). Since the solution is reliable only 
within the photosphere, the "foot" may not be a real feature. However, it is clear that 
when X, P > 0, the relative number of high column density lines of sight is diminished. For 
instance, the relative number of obscured AGNs per logarithm of column density reaches 
~ 1 with X = 0.08 at tt — 0.4, about half the number predicted with X = 0. The peak in 
the distribution at large column density could also be reduced if the inner edge were concave 
near the equatorial plane, i.e., if r in {z) were a decreasing function of \z\ close to the torus 
midplane. 
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CONCLUSIONS 



Previous work iKrolik! (120071 ) found the internal structure of radiation supported tori 
around AGN when they are heated only at their inner edge. In this paper, we have found 
generalized solutions to this problem by also taking into account two kinds of local heating 
sources, Compton scattering of hard X-rays and stars. Our results can be summarized in 
the following statements: 



1. For reasonable parameters, the local sources of heating can noticeably supplement the 
radiation pressure support to the tori. Local heating extends the matter distribution 
both radially and vertically. 

2. The effects of hard X-ray heating and stellar heating (at least in the Schmidt model) 
strongly resemble each other when their amplitudes are matched appropriately. 

3. Hydrostatic equilibrium in a torus with local heating demands a smaller range of 
(Lf Le)/tt than when there is none. It is ~ 0.1-0.2 for typical parameters. In addition, 
the equilibrium solutions require smaller orbital speed at the inner edge and a shallower 
gravitational potential inside. 

4. The local sources do little to change the predicted statistical distribution of column 
densities. 

5. In order to achieve hydrostatic equilibrium in the torus, the angular momentum has 
to be redistributed in both radial and vertical directions. 



Placing these formal results in context, we note that in a typical AGN we expect X < 
0.05 and P < 0.1. Because their effects add, we might expect that the effective amplitude of 
interior heating corresponds roughly to the X < 0.1 case. The data displayed in Figure [7^ 
would then imply a rather limited range of Q in which hydrostatic solutions can be found, 4 < 
Q <6- There are three possible conclusions that can follow from this inference: The first is 
that there are processes that automatically tune Q to lie in the permitted range. The second 
is that radiation pressure is so effective in supporting dusty gas against gravity that most 
tori are not in hydrostatic equilibrium. The third is that the numerous simplifications and 
approximations in our model have artificially narrowed the range of parameters permitting 
equilibria. We discuss them in this order. 

The central parameter combination controlling Q is (L/Le)/tt, for which the range 
corresponding to Q ~ 4-6 is ~ 0.1-0.15. The actual range of L/Le in AGN is not, at present 
well-known. Moreover, most AGN mass estimates done hitherto rely on the assumption that 
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radiation fo rces are unimportant compared to gravity. This is as true for those based on maser 
kinem atics (IGallimore et al.lll996l ; iGreenhill et al.l 119971 ; IVlemmings et all 120071 ; iKondratko 
20071 ) as for t hose based on broad emi s sion line widths and p hotoionization model-estimated 
lengthscales (IMcLure fc Dunlopll2004l ; iKollmeier et al.ll2006l ). It is the thrust of this paper, 
of course, that the dynamics of the molecular gas in which the masers are located might be 
substantially influenced by radiation forc es. In any event, the gen eral conclusion of these 
studies is that 0.1 < L/ L E < 1, although IMcLure fc Dunlopl (12004 ) might stretch the lower 
end of the range to ~ 0.03. Estimates of the characteristic Thomson depth t? are even harder 
to make, but the fact that a significant fraction of all type 2 AGN appear to have column 
densi ties for which tt > 1 (IRisaliti et al.1 Il999l ; lUeda et all 120071 ; iMartinez-Sansigre et al. 
20071 ) suggests that tt ~ 1, but with an unknown population dispersion, might be reasonable. 



Remarkably, given both the uncertainty in the observational estimates and the extremely 
simplified nature of the model presented here, the nominal range of Q predicted by our model 
agrees well with the range selected by observations if tt ~ 1. Correcting the observational 
range of L/Le for a possible systematic error due to the neglect of radiation forces would 
tend to move it to somewhat smaller values, which would, if anything, improve the match. 
However, even if the intrinsic b readth of the L/Le distribution is as little as ~ 4 (as advo- 
cated by IKollmeier et al.1 (120061 )). further tuning would be required if we take seriously the 
narrowness of our favored range for Q. One might imagine, for example, that tt adjusts 
in such a way as to put this ratio into the range permitting equilibrium (it is hard to see 
how, on the dynamical timescale of the torus, the Eddington ratio itself can be altered). 
While this might be possible, invoking such an effect begs the question of its mechanism: 
What causes the optical depth scale to change in precisely the way necessary to tune Q to a 
value permitting equilibrium? There might also be some partial loosening of the constraints 
due to variations in h. Smaller h at fixed Q would imply larger (L/L E )/r T , but also larger 
volumetric heating rates. 

On the other hand, failure of hydrostatic equilibrium due to excess radiation pressure 
raises other problems. If the radiation support is too large to be balanced, accretion fuel 
would be blown away, which might eventually lead to a reduction in L/L E and the possi- 
ble restoration, at least temporarily, of hydrostatic balance. However, there is a timescale 
mismatch problem: the mass- loss due to radiation occurs on a dynamical timescale (i.e., the 
orbital period), whereas inflow occurs much more slowly because it requires angular momen- 
tum transport. For this reason, readjustment of L/Le due to the expulsion of accretion fuel 
would be considerably slower than the fuel loss itself. For the same reason, resupply of the 
torus would be equally slow compared to the loss of torus material. Thus, a full-blown wind 
from the torus would lead to a long-lived state in which the nucleus continues to operate, 
but the torus is so depleted that there would be little obscuration. 
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The difficulties posed by lack of dynamical equilibrium in the torus can be made clearer 
by an explicit estimate of the associated mass-loss rate. Without knowing the dynamics of 
the outflow, we can roughly estimate the typical total mass outflow rate M ut by assuming 
that the escape speed is of the order of the local orbital speed and the outflow is isotropic. 
Under these assumptions, M ou t ~ M torus Vt ~ 80-Pr r\^ 2 M^ 2 M Q yr _1 , where the numerical 
value of Mtorus comes from our fiducial model, T\ = r» n /(l pc), and M 7 = M BH / (107M Q ). 
The accretion process of the black hole must then be very wasteful, as the accretion rate 
required to fuel the nucleus is only ~ 0.3(L/L E )M 7 (r]/0.1)~ 1 M Q yr^ 1 , where r\ is the usual 
radiative efficiency in rest-mass units. Indeed, it is even quite wasteful by the standards of 
the mass-loss rate that might be estimated on the basis of warm absorber column densities, 
~ IMj^Tw^c^sMq yr _1 , where r W)PC is the characteristic radius of the warm absorber 
outflow in parsecs and its column density is A^ in units of 10 23 H cm -2 . Moreover, as 
remarked in the previous paragraph, it is hard to see how such a large outflow rate could be 
maintained from the outside. Estimating the resupply rate in conventional a-model terms, 
we find Mi n ~ Sttt^ 2 M^ 2 ao,i(h /r) 2 M Q yr _1 , where h/r is the ratio of the scale height to 
the radius, and ao.i = a/0.1. 

Lastly, we consider the possibility that a more careful or complete calculation might lead 
to more precise agreement with the observed range of L/Le (and tt, when that is better 
measured). There are several improvements to our calculation that might well improve 
its quantitative accuracy: replacing a single averaged opacity with one that depends on 
frequency and substituting genuine transfer for the diffusion approximation are two that come 
immediately to mind. In addition, there is the intriguing possibility that a radiation-driven 
outflow (or possibly even a radiation-supported equilibrium) might be unstable to short- 
wavelength compr essive fluctuations, i.e., clumping. Such a res ult may also be attractive 
for other reasons (IKrolik fc Begelmanl Il988l ; iNenkova et al.l 120021 ) . If the clumping is strong 
enough to reduce the effective IR optical depth of the torus to < h, the radiation force 
would be diminished. This process might then be self- limiting, as the radiation dynamics 
responsible for clumping in the first place would likewise be weakened. Exploring these 
possibilities is, of course, an enterprise we must leave for future work. 



We thank Eliot Quataert, Norm Murray, and Phil Chang for conversations that helped 
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-19- 



REFERENCES 

Barthel, P. 1989,ApJ, 336,606 

Beckmann, V., Gehrels, N., Shrader, C.R., & Soldi, S. 2006, ApJ, 638, 642 
Chang, P., Quataert, E., & Murray, N. 2007, ApJ, 662, 94 



Davies, R.I. et al. 2007, ApJ, in presss( ]astro-ph/0704137| l) 

di Serego Alighieri, S., Cimatti, A., & Fosbury, R.A.E. 1994, ApJ, 431, 123 

Efstathiou, A., & Rowan-Robinson,M. 1995, MNRAS, 273, 649 

Gallimore, J.F., Baum, S.A., O'Dea, CP, Brinks, E. & Pedlar, A. 1996, ApJ, 462, 740 

Granato, G.L., Danese, L., & Franceschini, A. 1997, ApJ, 486, 147 

Greenhill, L.J., Moran, J.M. & Herrnstein, J.R. 1997, ApJLetts, 481, L23 

Hao, L. et al. 2005, AJ, 2005, 129, 1795 

Jaffe, W. et al. 2004, Nature, 429, 47 

Kennicutt, Robert C.,Jr. 1998, ApJ, 498, 541 

Kollmeier, J.A. et al. 2006, ApJ, 648, 128 

Kondratko, P.T. 2007, unpublished Harvard Ph.D. thesis 

Kdnigl, A., & Kartje, J.F. 1994, ApJ, 434, 446 

Krolik, J.H., & Begelman, M.C. 1988, ApJ, 329, 702 

Krolik, J.H. 1999, Active Galactic Nuclei: From the Central Black Hole to the Galactic 
Environment, Jeremiah P. Ostriker and David N. Spergel, eds., (Princeton: New 
Jersey) 245 

Krolik, J.H. 2007, ApJ, 661, 52 

Marconi, A. et al. 2004, MNRAS, 351, 169 

Martinez-Sansigre, A. et al. 2007, MNRAS, 379, L6 

McLure, R.J. & Dunlop, J.S. 2004, MNRAS, 352, 1390 

Nenkova, M., Ivezic, Z., & Elitzur, M. 2002, ApJL, 570, L9 



-20- 

Pier, E.A., & Krolik, J.H. 1992a, ApJ, 399, L23 

Pier, E.A., & Krolik, J.H. 1992b, ApJ, 401, 99 

Risaliti, D., Maiolino, R., & Salvati, M. 1999, ApJ, 522, 157 

Semenov, D. et al. 2003, A&A, 410, 611 

Thompson, T., Quataert, E., & Murray, N. 2005, ApJ, 630, 167 
Tozzi, R. et al. 2006, A&A, 451, 457 
Triester, E et al. 2004, ApJ, 616, 123 
Ueda, Y. et al. 2007, ApJ, 664, 79 

Vlemmings, W.H.T., Bignall, H.E., & Diamond, P.J. 2007, ApJ, 656, 198 
Zakamska, N. et al. 2005, AJ, 129, 1212 



This preprint was prepared with the AAS IATgX macros v5.2. 



-21 - 




Fig. 1.- 



Flow chart of the numerical approach described in section 13.41 
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Fig. 2 — Solution with j in = 0.5, a = 1.5, r* = 10, Q = 4, 7 = 1.43 and X = 0.02. Left: Radiation 
energy density. Right: Matter density. In both, the scale is logarithmic, and the thin white curves 
show the photospheres on the top of the torus. The white dashed line marks the radius outside 
of which we solve the combined hydrostatic and radiation diffusion equations; it is not a physical 
edge. 





Fig. 3.— Solution with j in = 0.5, a = 1.5, r* = 10, Q = 4, 7 = 1.43 and P = 2.5 x 10~ 2 . Left: 
Radiation energy density. Right: Matter density. In both, the scale is logarithmic, and the thin 
white curves show the photospheres on the top of the torus. The white dashed line marks the 
radius outside of which we solve the combined hydrostatic and radiation diffusion equations; it is 
not a physical edge. 
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Fig. 4. — Comparison between two solutions which are shown in Fig. [5] and El Left: Radiation 
energy density. Right: Matter density. Red solid lines: Stellar heating case; Blue dotted lines: 
X-ray heating; Green solid (dotted): Photosphere for stellar heating case (X-ray heating case). 
Contours are in logarithmic scale with separation of 0.2 on the left and 0.1 on the right. 
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Fig. 5. — Comparison between unperturbed solution (ji n = 0.5, a = 1.5, r* = 10, Q = 4.2 and 
7 = 1.5) and perturbed (same parameters as the unperturbed except Q = 4 and X = 0.06). Left: 
Radiation energy density. Right: Matter density. Red solid lines: Perturbed solution; Blue dotted 
lines: Unperturbed; Green dotted (solid): Unperturbed (perturbed) photosphere. Contours are in 
logarithmic scale with separation of 0.2 on the left and 0.1 on the right. 
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Fig. 6. — Distribution of j 2 for the solutions in FigHJ Red solid lines: Perturbed. Blue dotted: 
Unperturbed. Contours are scaled linearly with values 0.3 to 1.0 from left to right. 
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Fig. 7. — Solution distributions on (a) Q — 7, (b) ji n — 7 and (c) a — 7 planes with fixed other 
parameters, and (d) the relationship between P and the stellar luminosity. Typical values ji n = 0.5, 
a = 1.5, t* = 10 and Q = 4 are adopted if not specified. 
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Fig. 8. — Predicted column density distribution for solutions with ji n — 0.5, a = 1.5, r* = 10, 
tt = 0.5, and fixed L uv . Left: X-ray heating case (top); Stellar heating (bottom). Right: Enlarged 
sections of the left. See Table [JJ and [2] for descriptions of the lines and corresponding X/P, Q and 
7- 
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Table 1. Parameters in top two panels in Fig. [8] 



Line Color 


Line Style 


X 


Q 


7 


Black 


Solid 


0.0 


4.35 


1.68 


Green 


Dashed 


0.02 


4.26 


1.63 


Red 


Dash- Dot 


0.06 


4.10 


1.57 


Blue 


Long-Dashes 


0.08 


4.0 


1.53 



Table 2. Parameters in bottom two panels in Fig. [8] 



Line Color 


Line Style 


P 


Q 


7 


Black 


Solid 


0.0 


4.35 


1.68 


Green 


Dashed 


0.02 


4.25 


1.63 


Red 


Dash- Dot 


0.05 


4.10 


1.55 


Blue 


Long-Dashes 


0.07 


4.0 


1.51 



